Two day spinup from Levitus. Interestingly this case highlights the sensitivity of the various parameters to U velocities which are 0 at day 1 and just partially spun up at day 2. As constructed the Suppression parameter is structured by Zonal velocities and shows the most variablility.
Working from these pieces I tried a little experiment with the 1st day of data where I replaced the U_mean field with the surface ECCO velocity field, holding U_mean constant throughout the 100-2000m layer. The resulting Suppression field looks very much like the Bate's et al. and Kappa shows corresponding improvements in the tropics with the Ecco Velocities providing the zonal structure in the tropics. There is even a hint of the second zonal low just south of the equator.
\(\color{green} {Parameterizing \quad the\quad Eddy\quad Length\quad Scale} \)
NOTES:
\(K=u_{rms}∗{\Gamma * \color{red}{L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}\)
\(L_{eddy} = min (L_r, L_{req}, \require{cancel}\cancel{L_{Rh}})\)
Since the two length scales we are using in the parameterization of \(L_{eddy}\) depend on the Baroclinic wave speed \(c_r\) we will first check \(c_r\) against Chelton 1998.
\(\quad\quad\quad\quad\quad\quad\quad\) First Baroclinic Wave Speed \(c_r\) CESM
# <!-- collapse=True -->
### get Variable from file c_rossby is in cm/s convert to m/s
c_rossby = nc.variables['C_ROSSBY'][0,:,:]
#Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
targ_def25, radius_of_influence=500000, fill_value=None)
c_rossby_mps = c_rossby_nearest/100.
clevs=arange(1.,3.6,.2)
fig=MapZoneContour(targ_def25.lons,targ_def25.lats,c_rossby_mps,figsize=(14,9),
title="CESM First Baroclinic Wave Speed (m/s)",
levels=clevs,
llcrnrlon=25.,
fmt='%.1f',
imagename='Fig_1_CESM_First_Baroclinic_Wave_Speed.png')
#display(fig1)
#DPI = fig1.get_dpi()
#print "DPI:", DPI
#DefaultSize = fig1.get_size_inches()
#print "Default size in Inches", DefaultSize
#print "Which should result in a %i x %i Image"%(DPI*DefaultSize[0], DPI*DefaultSize[1])
# the default is 100dpi for savefig:
\(\quad\quad\quad\quad\quad\quad\) Chelton Baroclinic gravity wave phase speed
# <!-- collapse=True -->
c_r_chelton=Image(filename=basedir+'/steering/chelton_sfig1.jpg')
display(c_r_chelton)
\(\quad\quad\quad\) The barclinic wave speed looks reasonable so lets compare the Paramaterized Rossby Radius to Chelton
\(\quad\quad\quad\quad\quad\quad\) First Baroclinic Rossby Radius CESM \(\quad (L_{eddy})\)
# <!-- collapse=True -->
fcor = nc.variables['FCOR'][0,:,:]
fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
targ_def25, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
l_rossby_calc=c_rossby_nearest/absf
btp = nc.variables['BTP'][0,:,:]
btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
targ_def25, radius_of_influence=500000, fill_value=None)
l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)
### convert regridded radius from cm to km
cm2km=1.e-5
l_eddy_calc_km=l_eddy_calc*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig=MapZoneContour(targ_def25.lons,targ_def25.lats,l_eddy_calc_km,figsize=(14,9),
levels=clevs,
llcrnrlon=25.,
title="CESM First Baroclinic Rossby Radius (km) $L_{eddy} \sim min(l_r,l_{req})$",
fmt='%.1i',
imagename='Fig_2_CESM_First_Baroclinic_Rossby_Radius.png')
lrhines = nc.variables['L_RHINES'][0,:,:]
lrhines_nearest = pyresample.kd_tree.resample_nearest(orig_def, lrhines, \
targ_def25, radius_of_influence=500000, fill_value=None)
###clevs=np.logspace(1,15,40)
###norm=LogNorm()
###fig,ax,cbar =MapContourg(targ_def0,lrhines_nearest,
### addzonal=True,
### levels1=clevs,
### norm=norm,
### figsize=(14,5),
### title="$L_{rhines} (cmn^2)$ log scale")
###cbar.set_ticks(np.logspace(1,15,15))
if False:
l_eddy1=np.maximum(l_eddy_calc,lrhines_nearest)
l_eddy1_km=l_eddy1*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig2=MapZoneContour(targ_def25.lons,targ_def25.lats,l_eddy1_km,figsize=(14,9),
levels=clevs,
llcrnrlon=25.,
title="CESM First Baroclinic Rossby Radius (km) using rhines scale $L_{eddy} \sim max(l_rhines,min(l_r,l_{req}))$",
fmt='%.1i' )
# <!-- collapse=True -->
l_r_chelton=Image(filename=basedir+'/steering/chelton_fig2.jpg')
display(l_r_chelton)
# <!-- collapse=True -->
# recalculate these based on maps which start at 0 longitude
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
targ_def0, radius_of_influence=500000, fill_value=None)
c_rossby_mps = c_rossby_nearest/100.
fcor = nc.variables['FCOR'][0,:,:]
fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
targ_def0, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
l_rossby_calc=c_rossby_nearest/absf
btp = nc.variables['BTP'][0,:,:]
btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
targ_def0, radius_of_influence=500000, fill_value=None)
l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)
cm2km=1.e-5
l_eddy_calc_km=l_eddy_calc*cm2km
lrhines = nc.variables['L_RHINES'][0,:,:]
lrhines_nearest = pyresample.kd_tree.resample_nearest(orig_def, lrhines, \
targ_def0, radius_of_influence=500000, fill_value=None)
l_eddy1=np.maximum(l_eddy_calc,lrhines_nearest)
l_eddy1_km=l_eddy1*cm2km
\(\color{green} {Parameterizing \quad Zonal\quad Eddy\quad Phase\quad Speed}\quad (\color{red}{c})\)
NOTES:
\(K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - \color{red}{c}|^2 /u_{rms}^2 (z=0)}\)
\(\require{cancel}\cancel{c = - \beta * {L_r^2}}\) \(L_r\) too high at equator
\(c = - \beta * L_{eddy}^2\)
\(\quad\quad\quad\quad\quad\) Eddy Phase speed CESM
# <!-- collapse=True -->
if False:
#beta
clevs=arange(1.e-14,25.e-14,1e-14)
fig,ax,cbar =MapContourg(targ_def0,btp_nearest,
levels1=clevs,
# norm=norm,
addzonal=True,
figsize=(14,5),
title="Beta (1/cm*s)",
imagename='Fig_3_Beta.png')
l_eddy_calc_sq=l_eddy_calc*l_eddy_calc
if True:
clevs=np.logspace(5,8,40)
norm=LogNorm()
fig,ax,cbar =MapContourg(targ_def0,l_eddy_calc,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
title="$L_{eddy} (cm)$ log scale")
# imagename='Fig_4_L_eddy_SQ.png')
#cbar.ax.get_xaxis().get_major_formatter().labelOnlyBase = True
#cbar.ax.get_yaxis().set_major_formatter(plt.LogFormatter(10, labelOnlyBase=True))
#ax.yaxis.set_major_locator(ticker.MultipleLocator)
#ax.yaxis.set_major_formatter(LogFormatter())
#cbar.formatter.set_scientific(False)
#cbar.formatter.set_powerlimits((0, 0))
#cbar.update_ticks()
cbar.set_ticks(np.logspace(5,8,15))
if False:
clevs=np.logspace(1,15,40)
norm=LogNorm()
fig,ax,cbar =MapContourg(targ_def0,l_eddy_calc_sq,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
title="$L_{eddy}^2 (cm^2)$ log scale")
# imagename='Fig_4_L_eddy_SQ.png')
#cbar.ax.get_xaxis().get_major_formatter().labelOnlyBase = True
#cbar.ax.get_yaxis().set_major_formatter(plt.LogFormatter(10, labelOnlyBase=True))
#ax.yaxis.set_major_locator(ticker.MultipleLocator)
#ax.yaxis.set_major_formatter(LogFormatter())
#cbar.formatter.set_scientific(False)
#cbar.formatter.set_powerlimits((0, 0))
#cbar.update_ticks()
cbar.set_ticks(np.logspace(1,15,15))
fig.savefig('Fig_4_L_eddy_SQ.png',bbox_inches='tight',dpi=200)
c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s
clevs=arange(-10,20,1)
fig,ax,cbar = MapContourg(targ_def0,c,
addzonal=True,
levels1=clevs,
# norm=norm,
setover='darkred',
figsize=(14,5),
title="eddy phase speed c = beta*$l_{eddy}^2$ (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")
cbar.set_ticks(np.linspace(-10,20,7))
fig.savefig('Fig_5_eddy_phase_speed.png',bbox_inches='tight',dpi=200)
if False:
c_rhines=btp_nearest*l_eddy1*l_eddy1
clevs=arange(-10,20,1)
fig,ax,cbar = MapContourg(targ_def0,c_rhines,
addzonal=True,
levels1=clevs,
# norm=norm,
setover='darkred',
figsize=(14,5),
title="eddy phase speed including rhines scale c = beta*$l_{eddy}^2$ (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")
cbar.set_ticks(np.linspace(-10,20,7))
c=-1.*c
\(\quad \quad \quad \quad \)Compared to Hughes, the phase speed is much too high near the equator.
\(\quad \quad \quad \quad \)The first tuning mod is to limit the phase speed to 20cm/s.
\(\quad \quad \quad \quad \)Eddy Phase speed CESM \((c = max(- \beta * L_{eddy}^2,-20) )\)
# <!-- collapse=True -->
ceddy2 = nc.variables['C_EDDYLIM'][0,:,:]
ceddy2_nearest = pyresample.kd_tree.resample_nearest(orig_def, ceddy2, \
targ_def0, radius_of_influence=500000, fill_value=None)
###c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s
clevs=arange(-10,21,1)
fig,ax,cbar = MapContourg(targ_def0,-1*ceddy2_nearest,
addzonal=True,
levels1=clevs,
# norm=norm,
setover='darkred',
figsize=(14,5),
title="eddy phase speed c = beta*leddy**2 (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")
yticks=np.linspace(-10,20,7)
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
fig.savefig('Fig_6_eddy_phase_speed_linear.png',bbox_inches='tight',dpi=200)
\(\quad\quad\quad\quad\quad\quad\) Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09
# <!-- collapse=True -->
hughes_c=Image(filename=basedir+'/steering/Hughes_phase_speed.png')
display(hughes_c)
\(\color{green}{Zonal\quad Velocity \quad Term}\quad (\color{red}{u_{mean}})\)
NOTES:
\(K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |\color{red}{u_{mean}} - c|^2 /u_{rms}^2 (z=0))}\)
Here \(U_{mean}\) is just CESM surface velocity.
\(\quad\quad \quad \quad\quad \) Zonal Velocity CESM
# <!-- collapse=True -->
################# PLOT (u-c)^2 convert m^2/s^2 #######################
umean=nc.variables['U_MEAN'][0,:,:]
umean10_nearest = pyresample.kd_tree.resample_nearest(orig_def, umean,
targ_def0, radius_of_influence=500000, fill_value=None)
umean10_mps=umean10_nearest*1.e-2
clevs=np.linspace(-60,60,40)
fig,ax,cbar=MapContourg(targ_def0,umean10_nearest,
addzonal=True,
figsize=(14,5),
levels1=clevs,
# setover='darkred',
# setunder='darkblue',
title="U zonal velocity day 1 (cm/s)")
#fig.subplots_adjust(left=0.15, right=.9, bottom=0.1, top=0.9);
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.]);
fig.savefig('Fig_7_POP_UVEL.png',bbox_inches='tight',dpi=200)
utgrd_md1_grd0 = ncday1.variables['UTGRD_MAX'][0,:,:,:]
utgrd_md1_nearest = pyresample.kd_tree.resample_nearest(orig_def,utgrd_md1_grd0[0,:,:] ,
targ_def0, radius_of_influence=500000, fill_value=None)
clevs=np.linspace(-60,60,40)
fig,ax,cbar=MapContourg(targ_def0,utgrd_md1_nearest,
addzonal=True,
figsize=(14,5),
levels1=clevs,
# setover='darkred',
# setunder='darkblue',
title="U zonal velocity day 2 (cm/s)")
#fig.subplots_adjust(left=0.15, right=.9, bottom=0.1, top=0.9);
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.]);
if False:
# U averaged over 0-50m
clevs=np.linspace(-60,60,40)
fig,ax,cbar=MapContourg(targ_def0,u50_nearest,
addzonal=True,
figsize=(14,5),
levels1=clevs,
# setover='darkred',
# setunder='darkblue',
title="U 0-50m avg zonal velocity (cm/s)")
#fig.subplots_adjust(left=0.15, right=.9, bottom=0.1, top=0.9);
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.]);
# get 2D versions of the lat and lon variables add longitude start point here!
lon_ecco_orig, lat_ecco_orig = np.meshgrid(lon_ecco[:], lat_ecco[:])
orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_ecco_orig, lats=lat_ecco_orig)
ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_ecco_def, ecco_u, \
targ_def0, radius_of_influence=500000, fill_value=None)
ecco_u_nearest_masked=np.ma.masked_where(ecco_u_nearest<-10, ecco_u_nearest, copy=True)
ecco_u_cd0_grd0=pyresample.kd_tree.resample_nearest(targ_def0,ecco_u_nearest_masked*100,
orig_def, radius_of_influence=500000, fill_value=None)
print ecco_u_cd0_grd0.max(),ecco_u_cd0_grd0.min()
#print ecco_u_cyc_masked[:,0]
clevs=np.linspace(-60,60,40)
#ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_def, ecco_u, \
# targ_def, radius_of_influence=500000, fill_value=None)
fig,ax,cbar=MapContourg(targ_def0,ecco_u_nearest_masked*100,
addzonal=True,
figsize=(14,5),
levels1=clevs,
# extend='both',
title="ECCO U zonal velocity (cm/s)")
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.]);
fig.savefig('Fig_8_ECCO_UVEL.png',bbox_inches='tight',dpi=200)
48.6693 -67.6657
\(\color{green}{(U-c)\quad Term}\)
NOTES:
- \(c = max(- \beta * L_r^2,-20)\) )
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\) (U-c) CESM
# <!-- collapse=True -->
uminusc2 = nc.variables['U_MINUS_C'][0,0,:,:]
uminusc2_nearest = pyresample.kd_tree.resample_nearest(orig_def, uminusc2, \
targ_def0, radius_of_influence=500000, fill_value=None)
###uminusc2_mps=uminusc2_nearest*1.e-4
if False:
clevs=np.linspace(-75,75,50)
fig,ax,cbar=MapContourg(targ_def0,uminusc2_nearest,
addzonal=True,
# norm=norm,
levels1=clevs,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="$|U-c|$ (cm/s) (max = %i/min = %i) (c calculated using new eddy phase speed)"%(uminusc2_nearest.max(),uminusc2_nearest.min()))
yticks=linspace(-75,75,6)
zonalticks=np.linspace(-25,25,3)
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
pylab.xlim([-5,25])
pylab.xticks(zonalticks);
clevs=np.linspace(0,25,25)
fig,ax,cbar=MapContourg(targ_def0,uminusc2_nearest,
addzonal=True,
# norm=norm,
levels1=clevs,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="$|U-c|$ day 1 (cm/s) (max = %i/min = %i) "%(uminusc2_nearest.max(),uminusc2_nearest.min()))
yticks=linspace(0,25,6)
zonalticks=np.linspace(0,20,3)
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
pylab.xlim([0,25])
pylab.xticks(zonalticks);
fig.savefig('Fig_10_U_minus_c_new.png',bbox_inches='tight',dpi=200)
uminusc2d1 = ncday1.variables['U_MINUS_C'][0,0,:,:]
uminusc2d1_nearest = pyresample.kd_tree.resample_nearest(orig_def, uminusc2d1, \
targ_def0, radius_of_influence=500000, fill_value=None)
clevs=np.linspace(0,25,25)
fig,ax,cbar=MapContourg(targ_def0,uminusc2d1_nearest,
addzonal=True,
# norm=norm,
levels1=clevs,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="$|U-c|$ day 2 (cm/s) c limited to 20 $cm/s$ (max = %i/min = %i) "%(uminusc2d1_nearest.max(),uminusc2d1_nearest.min()))
yticks=linspace(0,25,6)
zonalticks=np.linspace(0,20,3)
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
pylab.xlim([0,25])
pylab.xticks(zonalticks);
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\) (U-c) Bates et al
\(\color{green}{(U-c)^2 \quad Term}\)
NOTES:
- \(c = max(- \beta * L_{eddy}^2,-20)\)
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (U-c)^2\) CESM
# <!-- collapse=True -->
uminusc2sq = nc.variables['U_MINUS_C_SQ'][0,0,:,:]
uminusc2sq_nearest = pyresample.kd_tree.resample_nearest(orig_def, uminusc2sq, \
targ_def0, radius_of_influence=500000, fill_value=None)
uminusc2sq_mps=uminusc2sq_nearest*1.e-4
if False:
clevs=np.logspace(-5,0,40)
# norm=matplotlib.colors.SymLogNorm(linthresh=1e-4, vmin=-1000, vmax=-.001)
norm=LogNorm()
fig,ax,cbar=MapContourg(targ_def0,uminusc2sq_mps,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
title=" ${|(U-c)|}^2$ using new $c_{eddy} {(m^2/s^2)}$ max = %g/min = %g "%(uminusc2sq_mps.max(),uminusc2sq_mps.min()))
# Customize y tick lables
yticks=[1e-5,3e-5,1e-4,3e-4,1e-3,3e-3,1e-2,3e-2,1e-1,3e-1,1]
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
pylab.xlim([0.,.1])
fig.savefig('Fig_12_U_minus_c_SQ_new_log.png',bbox_inches='tight',dpi=200)
clevs=np.linspace(0,.06,25)
fig,ax,cbar=MapContourg(targ_def0,uminusc2sq_mps,
addzonal=True,
levels1=clevs,
setover='darkred',
setunder='black',
figsize=(14,5),
title=" ${|(U-c)|}^2$ day 1 ${(m^2/s^2)}$ max = %g/min = %g "%(uminusc2sq_mps.max(),uminusc2sq_mps.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
fig.savefig('Fig_13_U_minus_c_SQ_new.png',bbox_inches='tight',dpi=200)
uminusc2_calc_sq_mps=uminusc2*uminusc2*1.e-4
if False:
clevs=np.linspace(0,.06,25)
fig,ax,cbar=MapContourg(targ_def0,uminusc2_calc_sq_mps,
addzonal=True,
levels1=clevs,
setover='darkred',
setunder='black',
figsize=(14,5),
title=" ${|(U-c)|}^2 {(m^2/s^2)}$ calc max = %g/min = %g "%(uminusc2_calc_sq_mps.max(),uminusc2_calc_sq_mps.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
uminusc2sqd1 = ncday1.variables['U_MINUS_C_SQ'][0,0,:,:]
uminusc2sqd1_nearest = pyresample.kd_tree.resample_nearest(orig_def, uminusc2sqd1, \
targ_def0, radius_of_influence=500000, fill_value=None)
uminusc2sqd1_mps=uminusc2sqd1_nearest*1.e-4
clevs=np.linspace(0,.06,25)
fig,ax,cbar=MapContourg(targ_def0,uminusc2sqd1_mps,
addzonal=True,
levels1=clevs,
setover='darkred',
setunder='black',
figsize=(14,5),
title=" ${|(U-c)|}^2$ day 2 ${(m^2/s^2)}$ max = %g/min = %g "%(uminusc2sqd1_mps.max(),uminusc2sqd1_mps.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
fig.savefig('Fig_13_U_minus_c_SQ_new.png',bbox_inches='tight',dpi=200)
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad (U-c)^2\) Bates et al
# <!-- collapse=True -->
bates_uminuscsq=Image(filename=basedir+'/steering/Bates_u_minus_c_sq.jpg')
display(bates_uminuscsq)
\(\color{green}{Parameterizing \quad Eady \quad Growth \quad Rate}\quad \color{red}{\sigma}\)
NOTES:
- The original derivation of \(\sigma_{vi}\) goes to 0 at the equator because of \(f\)
\(\quad\quad\quad\quad\quad\sigma_{vi} = {f \over \sqrt{R_i}}\)
- The new derivation of \(\sigma_{vi}\):
\(\quad\quad\quad\quad\quad R_i = {f^2 N^2 \over {\underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}}\quad\quad = \quad\quad {f^2N^2 \over m^4} \)
\(\quad\quad\quad\quad\quad \sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}\)
\(\quad\quad\quad\quad\quad \)Eady Growth Rate \((\sigma) \quad\) CESM
# <!-- collapse=True -->
n2md0 = nc.variables['N2_MAX'][0,:,:,:]
m4md0 = nc.variables['M4_MAX'][0,:,:,:]
zwmd0 = nc.variables['z_w'][:]
dzwmd0= nc.variables['dzw'][:]
kmtmd0 = nc.variables['KMT'][:,:]
sigma_avg1_cd0_grd0=np.where(np.less_equal(0,kmtmd0[:,:])&np.greater(m4md0[0,:,:], 0.),0.,0.)
depthcd0=np.where( np.less(kmtmd0[:,:],0)&np.greater(m4md0[0,:,:], 0.) ,0.,0.)
for i in range(60):
if (zwmd0[i] >= 1.e4 and zwmd0[i]<2.e5):
sigma_avg1_cd0_grd0=sigma_avg1_cd0_grd0+np.where(np.less_equal(i,kmtmd0[:,:])&np.greater(m4md0[i,:,:], 0.),(n2md0[i,:,:]/m4md0[i,:,:]*dzwmd0[i]),0.)
depthcd0=depthcd0+np.where(np.less_equal(i,kmtmd0[:,:])&np.greater(m4md0[i,:,:], 0.),dzwmd0[i],0.)
sigma_avg1_cd0_grd0=np.where(depthcd0>0.,sigma_avg1_cd0_grd0/depthcd0,0.)
sigma_avg1_cd0_grd0=np.where(sigma_avg1_cd0_grd0>0.,1./np.sqrt(sigma_avg1_cd0_grd0),0.)
sigma_avg1_cd0_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1_cd0_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
# <!-- collapse=True -->
n2md1 = ncday1.variables['N2_MAX'][0,:,:,:]
m4md1 = ncday1.variables['M4_MAX'][0,:,:,:]
zwmd1 = ncday1.variables['z_w'][:]
dzwmd1= ncday1.variables['dzw'][:]
kmtmd1 = ncday1.variables['KMT'][:,:]
sigma_avg1_cd1_grd0=np.where(np.less_equal(0,kmtmd1[:,:])&np.greater(m4md1[0,:,:], 0.),0.,0.)
depthcd1=np.where( np.less(0,kmtmd1[:,:])&np.greater(m4md1[0,:,:], 0.) ,0.,0.)
for i in range(60):
if (zwmd1[i] >= 1.e4 and zwmd1[i]<2.e5):
sigma_avg1_cd1_grd0=sigma_avg1_cd1_grd0+np.where(np.less(i,kmtmd1[:,:])&np.greater(m4md1[i,:,:], 0.),(n2md1[i,:,:]/m4md1[i,:,:]*dzwmd1[i]),0.)
depthcd1=depthcd1+np.where(np.less(i,kmtmd1[:,:])&np.greater(m4md1[i,:,:], 0.),dzwmd1[i],0.)
sigma_avg1_cd1_grd0=np.where(depthcd1>0.,sigma_avg1_cd1_grd0/depthcd1,0.)
sigma_avg1_cd1_grd0=np.where(sigma_avg1_cd1_grd0>0.,1./np.sqrt(sigma_avg1_cd1_grd0),0.)
sigma_avg1_cd1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1_cd1_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
# <!-- collapse=True -->
sigma_avg1_md0_grd0 = nc.variables['SIGMA_AVG1'][0,:,:]
sigma_avg1_md0_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1_md0_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
sigma_avg1_md1_grd0 = ncday1.variables['SIGMA_AVG1'][0,:,:]
sigma_avg1_md1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1_md1_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
clevs=np.logspace(-6.6,-2.3,80,base=2,endpoint=True)
norm=LogNorm(vmin=.01,vmax=.21)
fig,ax,cbar=MapContourg(targ_def0,sigma_avg1_md0_nearest*86400,
levels1=clevs,
# addzonal=True,
norm=norm,
figsize=(14,5),
setunder='darkblue',
setover='darkred',
title="Eady inverse time scale day 1 from model $\sqrt(M4/N2) (days^{-1}) \sim {{m^2} \over N}$ max = %g/min = %g/mean = %g "%(sigma_avg1_md0_nearest.max()*86400,sigma_avg1_md0_nearest.min()*86400,sigma_avg1_md0_nearest.mean()*86400))
yticks=[.0125,.025,.05,.1,.2]
cbar.set_ticks(yticks)
cbar.ax.set_yticklabels(["1/80","1/40","1/20","1/10","1/5"]);
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,0.025,.05,.1]);
clevs=np.logspace(-6.6,-2.3,80,base=2,endpoint=True)
norm=LogNorm(vmin=.01,vmax=.21)
fig,ax,cbar=MapContourg(targ_def0,sigma_avg1_md0_nearest*86400,
levels1=clevs,
norm=norm,
figsize=(14,5),
setunder='darkblue',
setover='darkred',
title="Eady inverse time scale day 1 from model $\sqrt(M4/N2) (days^{-1}) \sim {{m^2} \over N}$ max = %g/min = %g/mean = %g "%(sigma_avg1_md0_nearest.max()*86400,sigma_avg1_md0_nearest.min()*86400,sigma_avg1_md0_nearest.mean()*86400))
yticks=[.04,.08,.12,.16,.2]
cbar.set_ticks(yticks)
if False:
clevs=np.logspace(-6.6,-2.3,80,base=2,endpoint=True)
norm=LogNorm(vmin=.01,vmax=.21)
fig,ax,cbar=MapContourg(targ_def0,sigma_avg1_cd0_nearest*86400,
levels1=clevs,
norm=norm,
figsize=(14,5),
setunder='darkblue',
setover='darkred',
title="Eady inverse time scale day 1 calculated $\sqrt(M4/N2) (days^{-1}) \sim {{m^2} \over N}$ max = %g/min = %g/mean = %g "%(sigma_avg1_cd0_nearest.max()*86400,sigma_avg1_cd0_nearest.min()*86400,sigma_avg1_cd0_nearest.mean()*86400))
yticks=[.04,.08,.12,.16,.2]
cbar.set_ticks(yticks)
clevs=np.logspace(-6.6,-2.3,80,base=2,endpoint=True)
norm=LogNorm(vmin=.01,vmax=.21)
fig,ax,cbar=MapContourg(targ_def0,sigma_avg1_md1_nearest*86400,
levels1=clevs,
norm=norm,
figsize=(14,5),
setunder='darkblue',
setover='darkred',
title="Eady inverse time scale day 2 from model $\sqrt(M4/N2) (days^{-1}) \sim {{m^2} \over N}$ max = %g/min = %g/mean = %g "%(sigma_avg1_md1_nearest.max()*86400,sigma_avg1_md1_nearest.min()*86400,sigma_avg1_md1_nearest.mean()*86400))
yticks=[.04,.08,.12,.16,.2]
cbar.set_ticks(yticks)
if False:
clevs=np.logspace(-6.6,-2.3,80,base=2,endpoint=True)
norm=LogNorm(vmin=.01,vmax=.21)
fig,ax,cbar=MapContourg(targ_def0,sigma_avg1_cd1_nearest*86400,
levels1=clevs,
norm=norm,
figsize=(14,5),
setunder='darkblue',
setover='darkred',
title="Eady inverse time scale day 2 calculated $\sqrt(M4/N2) (days^{-1}) \sim {{m^2} \over N}$ max = %g/min = %g/mean = %g "%(sigma_avg1_cd1_nearest.max()*86400,sigma_avg1_cd1_nearest.min()*86400,sigma_avg1_cd1_nearest.mean()*86400))
yticks=[.04,.08,.12,.16,.2]
cbar.set_ticks(yticks)
<matplotlib.figure.Figure at 0x7f08a87cba90>
\(\quad\quad\quad\quad\quad \)Eady Growth Rate \((\sigma) \quad\) Bates et al
# <!-- collapse=True -->
#tulloch_eady=Image(filename=basedir+'/steering/tulloch_eady.png')
tulloch_eady=Image(filename=basedir+'/steering/TullochPhd_Eady.png')
display(tulloch_eady)
\(\color{green}{Parameterizing \quad RMS \quad eddy \quad velocity}\quad \color{red}{u_{rms}}\)
NOTES:
\(K=\color{red}{u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2(z=0)}\)
\(u_{rms} = alpha*{\sigma_{vi}}*L_{eddy}\)
- alpha (scaling constant) = 4
\(\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}\)
\(u_{rms}\) is limited to 5 cm/s $ max(u_{rms},5.)$
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Velocity \((u_{rms})\quad\) CESM
# <!-- collapse=True -->
alpha=4
cm2m=.01
#
urms_avg = nc.variables['URMS_AVG'][0,:,:]
urms_avg_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg,
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg_mps=urms_avg_nearest*cm2m
if False:
clevs=np.logspace(-3,1,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(targ_def0,urms_avg_mps,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="$u_{rms} ({m \over s})$ using old sigma log plot(alpha= %.1f,max = %.1f/min = %.2G) "%(alpha,urms_avg_mps.max(),urms_avg_mps.min()))
#customize colorbar ticks
yticks=np.logspace(-5,1,7)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.22])
pylab.xticks([0,.1,.2]);
clevs=np.linspace(0,.4,40)
fig,ax,cbar=MapContourg(targ_def0,urms_avg_mps,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="$u_{rms} ({m \over s})$ using old sigma (alpha=4,max = %.1f/min = %.1f) "%(urms_avg_mps.max(),urms_avg_mps.min()))
#customize colorbar ticks
yticks=[0,.1,.2,.3,.4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.22])
pylab.xticks([0,.1,.2]);
urms_avg1 = nc.variables['URMS_AVG1'][0,0,:,:]
urms_avg1_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1,
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg1_mps=urms_avg1_nearest*cm2m
if False:
clevs=np.logspace(-3,1,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(targ_def0,urms_avg1_mps,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="$u_{rms} ({m \over s})$ using new sigma log plot(alpha= %.1f,max = %.1f/min = %.2G) "%(alpha,urms_avg1_mps.max(),urms_avg1_mps.min()))
#customize colorbar ticks
yticks=np.logspace(-5,1,7)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.22])
pylab.xticks([0,.1,.2]);
fig.savefig('Fig_16_Urms_new_log.png',bbox_inches='tight',dpi=200)
clevs=np.linspace(0,.4,40)
fig,ax,cbar=MapContourg(targ_def0,urms_avg1_mps,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="$u_{rms}$ day 1 $({m \over s})$ using new sigma (alpha= %.1f,max = %.1f/min = %.2G) "%(alpha,urms_avg1_mps.max(),urms_avg1_mps.min()))
#customize colorbar ticks
yticks=[0,.1,.2,.3,.4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.22])
pylab.xticks([0,.1,.2]);
fig.savefig('Fig_17_Urms_new.png',bbox_inches='tight',dpi=200)
cm2m=.01
l_eddy_calc_m=l_eddy_calc*cm2m
urms_calc_mps=alpha*sigma_avg1_md1_nearest*l_eddy_calc_m
if False:
clevs=np.linspace(0,.4,40)
fig,ax,cbar=MapContourg(targ_def0,urms_calc_mps,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="$u_{rms} ({m \over s})$ calc value (alpha= %.1f,max = %.1f/min = %.1f) "%(alpha,urms_calc_mps.max(),urms_calc_mps.min()))
#customize colorbar ticks
yticks=[0,.1,.2,.3,.4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.22])
pylab.xticks([0,.1,.2]);
if False:
clevs=np.logspace(5,8,60)
norm=LogNorm()
# clevs=np.linspace(0,1e9,40)
fig,ax,cbar =MapContourg(targ_def0,l_eddy_calc,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
title="$L_{eddy} (cm)$ log scale")
# imagename='Fig_4_L_eddy_SQ.png')
#cbar.ax.get_xaxis().get_major_formatter().labelOnlyBase = True
#cbar.ax.get_yaxis().set_major_formatter(plt.LogFormatter(10, labelOnlyBase=True))
#ax.yaxis.set_major_locator(ticker.MultipleLocator)
#ax.yaxis.set_major_formatter(LogFormatter())
#cbar.formatter.set_scientific(False)
#cbar.formatter.set_powerlimits((0, 0))
#cbar.update_ticks()
cbar.set_ticks(np.logspace(1,15,15))
urms_avg1d1 = ncday1.variables['URMS_AVG1'][0,0,:,:]
urms_avg1d1_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1d1,
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg1d1_mps=urms_avg1d1_nearest*cm2m
clevs=np.linspace(0,.4,40)
fig,ax,cbar=MapContourg(targ_def0,urms_avg1d1_mps,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="$u_{rms}$ day 2 $({m \over s})$ (alpha= %.1f,max = %.1f/min = %.2G) "%(alpha,urms_avg1d1_mps.max(),urms_avg1d1_mps.min()))
#customize colorbar ticks
yticks=[0,.1,.2,.3,.4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.22])
pylab.xticks([0,.1,.2]);
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad u_{rms} \) Bates et al
# <!-- collapse=True -->
bates_urms=Image(filename=basedir+'/steering/bates2013-urms.png')
display(bates_urms)
\(\color{green}{u_{rms}^2 \quad Term} \quad\)
NOTES:
- \(K={u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /\color{red}{u_{rms}^2} (z=0)} \)
- \(u_{rms}^2\) is a surface value
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Velocity Squared \((u_{rms}^2)\quad\) CESM
#### <!-- collapse=True -->
urms_avg1_sq = nc.variables['URMS_AVG1_SQ'][0,0,:,:]
urms_avg1_sq_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1_sq,
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg1_sq_mps=urms_avg1_sq_nearest*1.e-4
###urms_avg1_sq_mps =urms_avg1_mps*urms_avg1_mps
urms_avg_sq = nc.variables['URMS_AVG_SQ'][0,:,:]
urms_avg_sq_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg_sq,
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg_sq_mps=urms_avg_sq_nearest*1.e-4
if False:
clevs=np.logspace(-5,1,42)
norm=LogNorm()
fig,ax,cbar=MapContourg(targ_def0,urms_avg_sq_mps,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="$u_{rms}^2 old sigma({m^2 \over s^2})$ log plot (alpha=%.1f,max = %.2G/min = %.2G) "%(alpha,urms_avg_sq_mps.max(),urms_avg_sq_mps.min()))
#customize colorbar ticks
yticks=np.logspace(-5,1,7)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
##clevs=np.linspace(0,.06,40)
clevs=np.linspace(.001,.06,40)
fig,ax,cbar=MapContourg(targ_def0,urms_avg_sq_mps,
addzonal=True,
levels1=clevs,
setover='darkred',
setunder='darkblue',
figsize=(14,5),
title="$u_{rms}^2 old sigma ({m^2 \over s^2})$ (alpha=%.1f,max = %.2G/min = %.2G) "%(alpha,urms_avg_sq_mps.max(),urms_avg_sq_mps.min()))
#customize colorbar ticks
yticks=[.001,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
clevs=np.logspace(-5,1,42)
norm=LogNorm()
fig,ax,cbar=MapContourg(targ_def0,urms_avg1_sq_mps,
addzonal=True,
levels1=clevs,
norm=norm,
figsize=(14,5),
setover='darkred',
setunder='darkblue',
title="$u_{rms}^2 new sigma({m^2 \over s^2})$ log plot (alpha=%.1f,max = %.2G/min = %.2G) "%(alpha,urms_avg1_sq_mps.max(),urms_avg1_sq_mps.min()))
#customize colorbar ticks
yticks=np.logspace(-5,1,7)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
fig.savefig('Fig_18_Urms_SQ_log.png',bbox_inches='tight',dpi=200)
##clevs=np.linspace(0,.06,40)
clevs=np.linspace(.001,.06,40)
fig,ax,cbar=MapContourg(targ_def0,urms_avg1_sq_mps,
addzonal=True,
levels1=clevs,
setover='darkred',
setunder='darkblue',
figsize=(14,5),
title="$u_{rms}^2$ day 1 $({m^2 \over s^2})$ (alpha=%.1f,max = %.2G/min = %.2G) "%(alpha,urms_avg1_sq_mps.max(),urms_avg1_sq_mps.min()))
#customize colorbar ticks
yticks=[.001,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
fig.savefig('Fig_19_Urms_SQ.png',bbox_inches='tight',dpi=200)
urms_calc_mps_sq=urms_calc_mps*urms_calc_mps
if False:
clevs=np.linspace(.001,.06,40)
fig,ax,cbar=MapContourg(targ_def0,urms_calc_mps_sq,
addzonal=True,
levels1=clevs,
setover='darkred',
setunder='darkblue',
figsize=(14,5),
title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%.1f,max = %.2G/min = %.2G) "%(alpha,urms_calc_mps_sq.max(),urms_calc_mps_sq.min()))
#customize colorbar ticks
yticks=[.001,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
urms_avg1_sqd1 = ncday1.variables['URMS_AVG1_SQ'][0,0,:,:]
urms_avg1_sqd1_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1_sqd1,
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg1_sqd1_mps=urms_avg1_sqd1_nearest*1.e-4
clevs=np.linspace(.001,.06,40)
fig,ax,cbar=MapContourg(targ_def0,urms_avg1_sq_mps,
addzonal=True,
levels1=clevs,
setover='darkred',
setunder='darkblue',
figsize=(14,5),
title="$u_{rms}^2$ day 2 $({m^2 \over s^2})$ (alpha=%.1f,max = %.2G/min = %.2G) "%(alpha,urms_avg1_sq_mps.max(),urms_avg1_sq_mps.min()))
#customize colorbar ticks
yticks=[.001,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06]);
fig.savefig('Fig_19_Urms_SQ.png',bbox_inches='tight',dpi=200)
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad u_{rms}^2 \) Bates et al
# <!-- collapse=True -->
urms_bates=Image(filename=basedir+'/steering/Bates_urms_sq.jpg')
display(urms_bates)
\(\color{green}{{(U-c)}^2 \over {u_{rms}^2}} \quad \color{green}{Term}\)
NOTES:
- \(K={u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * \color{red}{|u_{mean} - c|^2/{u_{rms}^2} (z=0)}} \)
- \(c = max(- \beta * L_{eddy}^2,-20)\)
- \(u_{rms} = max(u_{rms},5.)\)
- \(u_{rms}^2\) is a surface value
\(\quad\quad\quad\quad\quad\quad\) \({(U-c)}^2 \over {u_{rms}^2} \quad\) CESM
# <!-- collapse=True -->
umc2ourms2 = nc.variables['UMCSQ_OVR_URMSSQ'][0,0,:,:]
umc2ourms2_nearest = pyresample.kd_tree.resample_nearest(orig_def,umc2ourms2 ,
targ_def0, radius_of_influence=500000, fill_value=None)
umc2ourms2d1 = ncday1.variables['UMCSQ_OVR_URMSSQ'][0,0,:,:]
umc2ourms2d1_nearest = pyresample.kd_tree.resample_nearest(orig_def,umc2ourms2d1 ,
targ_def0, radius_of_influence=500000, fill_value=None)
if False:
clevs=np.logspace(-1,2,42)
norm=LogNorm()
fig,ax,cbar=MapContourg(targ_def0,umc2ourms2_nearest,
addzonal=True,
levels1=clevs,
figsize=(14,4),
norm=norm,
setover='darkred',
setunder='darkblue',
title="${(U-c)}^2 \over u_{rms}^2$ log plot")
#customize colorbar ticks
#yticks=clevs=np.logspace(-1,2,10)
yticks=[1e-1,2e-1,4e-1,1,2,4,1e1,2e1,4e1,10e1]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,8])
pylab.xticks([0,4,8]);
fig.savefig('Fig_20_U_minus_c_ovr_Urms_SQ_log.png',bbox_inches='tight',dpi=200)
clevs=arange(0,5,.1)
fig,ax,cbar=MapContourg(targ_def0,umc2ourms2_nearest,
addzonal=True,
levels1=clevs,
figsize=(14,4),
setover='darkred',
setunder='darkblue',
title="${(U-c)}^2 \over u_{rms}^2$ day 1")
#customize colorbar ticks
yticks=[0,2,4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,8])
pylab.xticks([0,4,8]);
clevs=arange(0,5,.1)
fig,ax,cbar=MapContourg(targ_def0,umc2ourms2d1_nearest,
addzonal=True,
levels1=clevs,
figsize=(14,4),
setover='darkred',
setunder='darkblue',
title="${(U-c)}^2 \over u_{rms}^2$ day 2")
#customize colorbar ticks
yticks=[0,2,4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,8])
pylab.xticks([0,4,8]);
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad {|u-c|^2 \over u_{rms}^2} \) Bates et al
# <!-- collapse=True -->
bates_uminusc_sq_ovr_urms_sq =Image(filename=basedir+'/steering/Bates_suppress.jpg')
display(bates_uminusc_sq_ovr_urms_sq )
\(\color{green}{\quad Suppression \quad factor = {1 \over (1 + b1 * |\bar u - c|^2 /u_{rms (z=0)}^2 )}}\)
NOTES:
- \(c = max(- \beta * L_{eddy}^2,-20)\)
- \(u_{rms} = max(u_{rms},5.)\)
- \(u_{rms}^2\) is a surface value
- \(b1\) (scaling constant) = 4.
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor CESM
# <!-- collapse=True -->
supp = nc.variables['SUPP'][0,0,:,:]
supp_nearest = pyresample.kd_tree.resample_nearest(orig_def,supp ,
targ_def0, radius_of_influence=500000, fill_value=None)
clevs=np.linspace(0,1,30)
norm=matplotlib.colors.Normalize()
fig,ax,cbar=MapContourg(targ_def0,supp_nearest,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Suppression (limiting $u_{rms}$ to 5cm/s)= ${1 \over (1 + b1 * | u - c|^2 /(u_{rms (z=0)}^2)}$ day 1")
#customize colorbar ticks
yticks=np.linspace(0, 1, 6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.8])
pylab.xticks([0,.25,.5,.75]);
suppd1 = ncday1.variables['SUPP'][0,0,:,:]
suppd1_nearest = pyresample.kd_tree.resample_nearest(orig_def,suppd1 ,
targ_def0, radius_of_influence=500000, fill_value=None)
clevs=np.linspace(0,1,30)
norm=matplotlib.colors.Normalize()
fig,ax,cbar=MapContourg(targ_def0,suppd1_nearest,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Suppression (limiting $u_{rms}$ to 5cm/s)= ${1 \over (1 + b1 * | u - c|^2 /(u_{rms (z=0)}^2)}$ day 2 ")
#customize colorbar ticks
yticks=np.linspace(0, 1, 6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.8])
pylab.xticks([0,.25,.5,.75]);
# <!-- collapse=True -->
kappa_min=100.e4
kappa_max=10000.e4
b1=4.
gamma=1.75
alpha=4
cm2m=.01
kverttmp= ncregrid.variables['KVERT'][0,:,:,:]
kvert3d_nearest_raw=np.ma.array(np.empty(kverttmp.shape))
kvert3d_nearest_scalen2=np.ma.array(np.empty(kverttmp.shape))
z_t= ncregrid.variables['z_t'][:]
latregrid=ncregrid.variables['lat'][:]
lonregrid=ncregrid.variables['lon'][:]
lon=lonregrid
urms_avg1_md0_grd0 = nc.variables['URMS_AVG1'][0,0,:,:]
urms_avg1_md0_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1_md0_grd0,
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg1_sq_md0_grd0 = nc.variables['URMS_AVG1_SQ'][0,0,:,:]
urms_avg1_sq_md0_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1_sq_md0_grd0,
targ_def0, radius_of_influence=500000, fill_value=None)
utgrd_md0_grd0 = nc.variables['UTGRD_MAX'][0,:,:,:]
utgrd = ncregrid.variables['UTGRD_MAX'][0,:,:,:]
#utgrd = pyresample.kd_tree.resample_nearest(orig_def, utgrd_max,
# targ_def0, radius_of_influence=500000, fill_value=None)
n2norm_md0_grd0=nc.variables['BFREQ_SQ_NORM'][0,:,:,:]
lrhines_cd0_grd0=sigma_avg1_cd0_grd0/btp
absf_grd0=np.abs(fcor)
l_rossby_calc_grd0=c_rossby/absf_grd0
l_rossbyeq_calc_grd0=np.sqrt(c_rossby/(2*btp))
l_eddy_calc_grd0=np.minimum(l_rossby_calc_grd0,l_rossbyeq_calc_grd0)
l_eddy_calc_grd0_sq=l_eddy_calc_grd0*l_eddy_calc_grd0
urms_avg1_cd0_grd0=alpha*sigma_avg1_cd0_grd0*l_eddy_calc_grd0
urms_avg1_cd0_grd0=np.maximum(urms_avg1_cd0_grd0,5.)
urms_avg1_cd0_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1_cd0_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
urms_avg1_sq_cd0_grd0 = urms_avg1_cd0_grd0**2
urms_avg1_sq_cd0_nearest = pyresample.kd_tree.resample_nearest(orig_def, urms_avg1_sq_cd0_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
c_eddy_calc_grd0=-1.*btp*l_eddy_calc_grd0_sq
c_eddy_lim=np.maximum(c_eddy_calc_grd0,-20.)
kvert_cd0_grd0=np.where(np.greater(m4md0[:,:,:], 0.) ,0.,0.)
for i in range(60):
if i==0:
n2scale_md0_grd0=1.
else:
n2scale_md0_grd0=n2norm_md0_grd0[i,:,:]
urms_avg1_k_cd0_grd0=urms_avg1_cd0_grd0*n2scale_md0_grd0
urms_avg1_k_sq_cd0_grd0=urms_avg1_k_cd0_grd0*urms_avg1_k_cd0_grd0
# uminusc_cd0_grd0=np.where(np.less(i,kmtmd0[:,:])&np.greater(urms_avg1_sq_cd0_grd0[:,:], 0.),np.absolute(utgrd_md0_grd0[i,:,:]-c_eddy_lim[:,:]),0.)
uminusc_cd0_grd0=np.where(np.less(i,kmtmd0[:,:])&np.greater(urms_avg1_sq_cd0_grd0[:,:], 0.),np.absolute(ecco_u_cd0_grd0[:,:]-c_eddy_lim[:,:]),0.)
uminusc_cd0_grd0_masked=np.ma.masked_where(uminusc_cd0_grd0>100,uminusc_cd0_grd0, copy=True)
uminusc_cd0_grd0_masked=np.where(np.greater(uminusc_cd0_grd0,100),0.,uminusc_cd0_grd0)
umcsq_cd0_grd0=uminusc_cd0_grd0_masked*uminusc_cd0_grd0_masked
umcsq_ov_urmssq_cd0_grd0=np.where(np.less(i,kmtmd0[:,:])&np.greater(urms_avg1_sq_cd0_grd0[:,:], 0.),umcsq_cd0_grd0[:,:]/urms_avg1_sq_cd0_grd0[:,:],0.)
supp_cd0_grd0=np.where(np.less(i,kmtmd0[:,:])&np.greater(urms_avg1_sq_cd0_grd0[:,:], 0.),1./(1.+b1*umcsq_ov_urmssq_cd0_grd0[:,:]),0.)
lmix_cd0_grd0=np.where(np.less(i,kmtmd0[:,:])&np.greater(urms_avg1_sq_cd0_grd0[:,:], 0.),gamma*l_eddy_calc_grd0[:,:]*supp_cd0_grd0[:,:],0.)
kvert_cd0_grd0[i,:,:]=np.where(np.less(i,kmtmd0[:,:])&np.greater(urms_avg1_sq_cd0_grd0[:,:], 0.),urms_avg1_k_cd0_grd0[:,:]*lmix_cd0_grd0[:,:],0.)
kvert_cd0_grd0[i,:,:]=np.minimum(kvert_cd0_grd0[i,:,:],kappa_max)
kvert_cd0_grd0[i,:,:]=np.maximum(kvert_cd0_grd0[i,:,:],kappa_min)
kvertsrf_lim=kvert_cd0_grd0[0,:,:]
kvert3d_nearest_raw[i,:,:] = pyresample.kd_tree.resample_nearest(orig_def, kvert_cd0_grd0[i,:,:], \
targ_def0, radius_of_influence=500000, fill_value=None)
kvert3d_nearest_scalen2[i,:,:] =pyresample.kd_tree.resample_nearest(orig_def, kvertsrf_lim*n2scale_md0_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
if i==0:
supp_cd0_nearest = pyresample.kd_tree.resample_nearest(orig_def, supp_cd0_grd0, \
targ_def0, radius_of_influence=500000, fill_value=None)
clevs=np.linspace(0,1,30)
norm=matplotlib.colors.Normalize()
fig,ax,cbar=MapContourg(targ_def0,supp_cd0_nearest,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Suppression using Seasonal Ecco Velocities along with other day 0 parameters")
#customize colorbar ticks
yticks=np.linspace(0, 1, 6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.8])
pylab.xticks([0,.25,.5,.75]);
# kvertlim_cd0_grd0=np.minimum(kvert_cd0_grd0,kappa_max)
# kvertlim_cd0_grd0=np.maximum(kvertlim_cd0_grd0,kappa_min)
# kvert_cd0_nearest = pyresample.kd_tree.resample_nearest(orig_def, kvertlim_cd0_grd0[i,:,:], \
# targ_def0, radius_of_influence=500000, fill_value=None)
# kvert_cd0_nearest_mps= kvert_cd0_nearest*1.e-4
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor Bates
# <!-- collapse=True -->
bates_suppression =Image(filename=basedir+'/steering/bates2013_suppression.png')
display(bates_suppression )
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor (Zonal x depth) CESM
clevs=np.linspace(0,1,11)
z_t= ncregrid.variables['z_t'][:]
latregrid=ncregrid.variables['lat'][:]
lonregrid=ncregrid.variables['lon'][:]
lon=lonregrid
suppregrid= ncregrid.variables['SUPP'][0,:,:,:]
suppregridz=suppregrid.mean(axis=2)
clevs=np.linspace(0,1,100,endpoint=True)
yticks=[0,.2,.4,.6,.8,1.]
fig,ax,cbar=plotzonalocn(latregrid, z_t, suppregridz,
levels1=clevs,
levels2=clevs,
# equalcontourscale=True,
figsize=(10.5,6),
title="Suppression factor day 1 (custom contours max = %.1f/min = %.1f) "%(suppregridz.max(),suppregridz.min()))
# Add axis labels and plot title
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
# Customize y tick lables
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('fraction')
suppregridd1= ncregridd1.variables['SUPP'][0,:,:,:]
suppregridzd1=suppregridd1.mean(axis=2)
clevs=np.linspace(0,1,100,endpoint=True)
yticks=[0,.2,.4,.6,.8,1.]
fig,ax,cbar=plotzonalocn(latregrid, z_t, suppregridzd1,
levels1=clevs,
levels2=clevs,
# equalcontourscale=True,
figsize=(10.5,6),
title="Suppression factor day 2 (custom contours max = %.1f/min = %.1f) "%(suppregridzd1.max(),suppregridzd1.min()))
# Add axis labels and plot title
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
# Customize y tick lables
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('fraction')
\(\quad\quad\quad\quad\quad\quad\quad\) Suppression Factor (Zonal x depth) Bates
# <!-- collapse=True -->
urms_bates7=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-007.jpg')
display(urms_bates7 )
\(\color{green} {\quad Parameterizing \quad The \quad Mixing \quad Term \quad }\color{red}{L_{mix}}\)
NOTES:
\(\color{red}{L_{mix}} = \Gamma * L_{eddy} * Suppression\)
\(Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}\)
\(\color{red}{L_{mix}} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}\)
\(\Gamma = 1.75\) (Tuning mod: the original Gamma of .35 produced a Kappa with the correct structure but too weak.
\(\quad\quad\quad\quad\quad\quad\quad\) LMIX CESM
lmix2 = nc.variables['LMIX'][0,0,:,:]
lmix2_nearest = pyresample.kd_tree.resample_nearest(orig_def,lmix2*cm2m ,
targ_def0, radius_of_influence=500000, fill_value=None)
norm=LogNorm()
clevs=np.logspace(0,5,40)
fig,ax,cbar=MapContourg(targ_def0,lmix2_nearest,
addzonal=True,
norm=norm,
levels1=clevs,
figsize=(14,5),
title="$L_{mix} = {\Gamma * L_{eddy} * Suppression} $ day 1")
cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,3e4,1e5])
#customize zonal ticks
##pylab.xlim([0,120e3])
##pylab.xticks([0,2500,5000,7500,10000]);
pylab.xlim([0,120e3])
##pylab.xticks([0,2500,5000,7500,10000]);
pylab.xticks([0,40000,80000,120000]);
lmix2d1 = ncday1.variables['LMIX'][0,0,:,:]
lmix2d1_nearest = pyresample.kd_tree.resample_nearest(orig_def,lmix2d1*cm2m ,
targ_def0, radius_of_influence=500000, fill_value=None)
norm=LogNorm()
clevs=np.logspace(0,5,40)
fig,ax,cbar=MapContourg(targ_def0,lmix2d1_nearest,
addzonal=True,
norm=norm,
levels1=clevs,
figsize=(14,5),
title="$L_{mix} = {\Gamma * L_{eddy} * Suppression} $ day 2")
cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,3e4,1e5])
#customize zonal ticks
##pylab.xlim([0,120e3])
##pylab.xticks([0,2500,5000,7500,10000]);
pylab.xlim([0,120e3])
##pylab.xticks([0,2500,5000,7500,10000]);
pylab.xticks([0,40000,80000,120000]);
\(\color{green}{\quad Parameterizing \quad Eddy \quad Diffusivity \quad} (\color{red}{K})\)
NOTES:
- \(\color{red}{K}=u_{rms}∗L_{mix}\)
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Diffusivity (K) CESM
# <!-- collapse=True -->
###gamma=.35
###############GAMMA bumped up for tuning #########################
gamma=1.75
b1=4.
###regrid2rect(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.017.pop.h.nday1.0249.nc',basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.017.regrid.nc',clobber=True)
###ncregrid = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.017.regrid.nc')
kvertsrf= ncregrid.variables['KVERT'][0,0,:,:]*1.e-4
clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(targ_def0,kvertsrf,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Eddy diffusivity day 1 $K$ (max = %.1f/min = %.1f) "%(kvertsrf.max(),kvertsrf.min()))
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3]);
kvertsrfd1= ncregridd1.variables['KVERT'][0,0,:,:]*1.e-4
clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(targ_def0,kvertsrfd1,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Eddy diffusivity day 2 $K$ (max = %.1f/min = %.1f) "%(kvertsrfd1.max(),kvertsrfd1.min()))
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3]);
clevs=np.linspace(0,1.e4,40)
kvert3d_nearest_lim_mps=kvert3d_nearest_raw[0,:,:]/1.e4
fig,ax,cbar=MapContourg(targ_def0,kvert3d_nearest_lim_mps,
addzonal=True,
levels1=clevs,
figsize=(14,5),
title="Eddy diffusivity using Ecco Velocities (max = %.1f/min = %.1f) "%(kvert3d_nearest_lim_mps.max(),kvert3d_nearest_lim_mps.min()))
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3]);
\(\quad\quad\quad\quad\quad\quad\quad\) Eddy Diffusivity (K) Bates
# <!-- collapse=True -->
bates_k =Image(filename=basedir+'/steering/bates2013_K.png')
display(bates_k )
\(\quad\quad\quad\quad\quad\quad\quad\) (Zonal x Depth) CESM
\(\quad\quad\quad\quad\quad\quad\quad\quad\quad\) K limited to (100 < K < 10000)
# <!-- collapse=True -->
###ncregrid = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.017.regrid.nc')
z_t= ncregrid.variables['z_t'][:]
latregrid=ncregrid.variables['lat'][:]
lonregrid=ncregrid.variables['lon'][:]
lon=lonregrid
kvert= ncregrid.variables['KVERT'][0,:,:,:]
kvertz=kvert.mean(axis=2)/1.e4
if False:
clevs=np.arange(0,100,20)
clevs=np.append(clevs,np.arange(100,1000,100))
clevs=np.append(clevs,np.arange(1000,4000,1000))
fig,ax,cbar=plotzonalocn(latregrid, z_t, kvertz,
levels1=clevs,
levels2=clevs,
equalcontourscale=True,
figsize=(10.5,6),
title="Eddy diffusivity day 1 $K$ (custom contours max = %.1f/min = %.1f) "%(kvertz.max(),kvertz.min()))
# Add axis labels and plot title
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
# Customize y tick lables
yticks=clevs
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('m**2/s')
clevs=np.linspace(0,2500,100,endpoint=True)
yticks=[0,500,1000,1500,2000,2500]
fig,ax,cbar=plotzonalocn(latregrid, z_t, kvertz,
levels1=clevs,
levels2=clevs,
# equalcontourscale=True,
figsize=(10.5,6),
title="Eddy diffusivity day 1 $K$ (Linear Scale max = %.1f/min = %.1f) "%(kvertz.max(),kvertz.min()))
# Add axis labels and plot title
levs=[500,1000,1500]
cs=ax.contour(latregrid, z_t, kvertz,levels=levs,colors='w',linewidths=2)
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
plt.clabel(cs, inline=1, fontsize=15)
# Customize y tick lables
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('m**2/s')
kvertd1= ncregridd1.variables['KVERT'][0,:,:,:]
kvertzd1=kvertd1.mean(axis=2)/1.e4
if False:
clevs=np.arange(0,100,20)
clevs=np.append(clevs,np.arange(100,1000,100))
clevs=np.append(clevs,np.arange(1000,4000,1000))
fig,ax,cbar=plotzonalocn(latregrid, z_t, kvertzd1,
levels1=clevs,
levels2=clevs,
equalcontourscale=True,
figsize=(10.5,6),
title="Eddy diffusivity day 2 $K$ (custom contours max = %.1f/min = %.1f) "%(kvertzd1.max(),kvertzd1.min()))
# Add axis labels and plot title
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
# Customize y tick lables
yticks=clevs
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('m**2/s')
clevs=np.linspace(0,2500,100,endpoint=True)
yticks=[0,500,1000,1500,2000,2500]
fig,ax,cbar=plotzonalocn(latregrid, z_t, kvertzd1,
levels1=clevs,
levels2=clevs,
# equalcontourscale=True,
figsize=(10.5,6),
title="Eddy diffusivity day 2 $K$ (Linear Scale max = %.1f/min = %.1f) "%(kvertzd1.max(),kvertzd1.min()))
# Add axis labels and plot title
levs=[500,1000,1500]
cs=ax.contour(latregrid, z_t, kvertzd1,levels=levs,colors='w',linewidths=2)
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
plt.clabel(cs, inline=1, fontsize=15)
# Customize y tick lables
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('m**2/s')
kvertz=kvert3d_nearest_raw.mean(axis=2)/1.e4
if False:
clevs=np.logspace(0,5,100,endpoint=True)
yticks=np.logspace(0,5,7,endpoint=True)
norm=LogNorm()
fig,ax,cbar=plotzonalocn(latregrid, z_t, kvertz,
levels1=clevs,
levels2=clevs,
norm=norm,
# equalcontourscale=True,
figsize=(10.5,6),
title="Eddy diffusivity from Ecco Velocities $K$ (Linear Scale max = %.1f/min = %.1f) "%(kvertz.max(),kvertz.min()))
# Add axis labels and plot title
levs=[500,1000,1500]
cs=ax.contour(latregrid, z_t, kvertz,levels=levs,colors='w',linewidths=2)
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
plt.clabel(cs, inline=1, fontsize=15)
# Customize y tick lables
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('m**2/s')
if False:
clevs=np.linspace(0,2500,100,endpoint=True)
yticks=[0,500,1000,1500,2000,2500]
fig,ax,cbar=plotzonalocn(latregrid, z_t, kvertz,
levels1=clevs,
levels2=clevs,
# equalcontourscale=True,
figsize=(10.5,6),
title="Eddy diffusivity using Ecco Velocities using model calc (max = %.1f/min = %.1f) "%(kvertz.max(),kvertz.min()))
# Add axis labels and plot title
levs=[500,1000,1500]
cs=ax.contour(latregrid, z_t, kvertz,levels=levs,colors='w',linewidths=2)
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
plt.clabel(cs, inline=1, fontsize=15)
# Customize y tick lables
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('m**2/s')
kvertzlim1=kvert3d_nearest_scalen2.mean(axis=2)/1.e4
clevs=np.linspace(0,2500,100,endpoint=True)
yticks=[0,500,1000,1500,2000,2500]
fig,ax,cbar=plotzonalocn(latregrid, z_t, kvertzlim1,
levels1=clevs,
levels2=clevs,
# equalcontourscale=True,
figsize=(10.5,6),
title="Eddy diffusivity using Ecco velocities scaled from surface value limited 100:10000 (max = %.1f/min = %.1f) "%(kvertzlim1.max(),kvertzlim1.min()))
# Add axis labels and plot title
levs=[500,1000,1500]
cs=ax.contour(latregrid, z_t, kvertzlim1,levels=levs,colors='w',linewidths=2)
plt.xlabel('Latitude')
plt.ylabel('Depth (cm)')
plt.clabel(cs, inline=1, fontsize=15)
# Customize y tick lables
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('m**2/s')
Here is the Bate's version
urms_bates8=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-008.jpg')
display(urms_bates8)
\(\quad\quad\quad\quad\quad\quad\quad\) Zonal average of the N2 normalized scaling CESM
clevs=np.linspace(0,1,11)
###ncregrid = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.015.regrid.nc')
n2norm=ncregrid.variables['BFREQ_SQ_NORM'][0,:,:,:]
n2normz=n2norm.mean(axis=2)
fig,ax,cbar=plotzonalocn(latregrid, z_t, n2normz,
levels1=clevs,
levels2=clevs,
equalcontourscale=True,
figsize=(14,5),
title="N2 Scaling day 1(max = %.1f/min = %.1f) "%(n2normz.max(),n2normz.min()))
# Add axis labels and plot title
plt.xlabel('Latitude')
plt.ylabel('Scaling factor')
# Customize y tick lables
yticks=clevs
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('fraction')
clevs=np.linspace(0,1,11)
###ncregrid = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.015.regrid.nc')
n2normd1=ncregridd1.variables['BFREQ_SQ_NORM'][0,:,:,:]
n2normzd1=n2normd1.mean(axis=2)
fig,ax,cbar=plotzonalocn(latregrid, z_t, n2normzd1,
levels1=clevs,
levels2=clevs,
equalcontourscale=True,
figsize=(14,5),
title="N2 Scaling day 2(max = %.1f/min = %.1f) "%(n2normzd1.max(),n2normzd1.min()))
# Add axis labels and plot title
plt.xlabel('Latitude')
plt.ylabel('Scaling factor')
# Customize y tick lables
yticks=clevs
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
cbar.set_label('fraction')
fig.savefig('Fig_31_N2_NORM_zonal_depth.png',bbox_inches='tight',dpi=200)
# <!-- collapse=True -->
from IPython import utils
from IPython.core.display import HTML
import os
def css_styling():
"""Load default custom.css file from ipython profile"""
base = utils.path.get_ipython_dir()
styles = "<style>\n%s\n</style>" % (open(os.path.join(base,'profile_default/static/custom/custom.css'),'r').read())
return HTML(styles)
css_styling()
The alternate derivation of \(\sigma_{vi}\):
\(R_i = {N^2\over{(\frac{\partial u}{\partial z})^2+(\frac{\partial v}{\partial z})^2}}\)
\(N^2={-g \over \rho_0 }\frac{\partial \rho}{\partial z}\)
After hydrostatic and geostrophic approximations
\(f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} \)
so
\(\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}}\)
\(R_i = {f^2 N^2 \over {\underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}}\quad\quad = \quad\quad {f^2N^2 \over m^4} \)
\(\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}\)
\(RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}\)
\(RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}\)
\(RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}\)
\(\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-g\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz\)
Note: missing \(f^2\) which will be cancelled when forming \(\sigma_{vi}\)
\(\quad\quad\) so \(\cdots\) this is not \(R_i\)
Implementation notes
Numerator : Top \(= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)\)
Denominator :
\(\begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align}\)
\(\begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 / DYT(i,j)^2 \\ \end{align}\)
\(\begin{align} work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k) \end{align}\)
Notes:
1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2